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Three-dimensional  (3D)  multi-physics  models  of  co-,  counter-  and  cross-flow  planar  solid  oxide  fuel  cell 
(SOFC)  stack  units  are  described.  The  models  consider  electronic  conduction  in  the  electrodes,  ionic  con¬ 
duction  in  the  electrolyte,  mass  transport  in  the  porous  electrodes  and  electrochemical  reactions  on  the 
three  phase  boundaries.  Based  on  the  analysis  of  the  ionic  conducting  equation  for  the  thin  electrolyte 
layer,  a  mathematically  equivalent  method  is  proposed  to  scale  the  electrolyte  thickness  with  the  cor¬ 
responding  change  in  the  ionic  conductivity  to  moderate  the  thin  film  effect  in  the  meshing  step  and 
decrease  the  total  number  of  degrees  of  freedom  in  the  3D  numerical  models.  Examples  of  applications 
are  given  with  typical  physical  fields  illustrated  and  the  characteristic  features  discussed  for  co-,  counter- 
and  cross-flow  designs.  The  3D  models  are  also  used  to  optimize  the  rib  widths  in  SOFC  stacks  as  a  function 
of  interconnect-electrode  contact  resistance. 

©  2009  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

Anode-supported  planar  solid  oxide  fuel  cell  (SOFC)  has  received 
much  attention  in  recent  years  due  to  its  satisfactory  power  density 
at  the  intermediate  temperature  [1-4].  The  reduced  operation  tem¬ 
perature  of  planar  SOFC  also  enables  cheaper  materials  to  be  used 
as  interconnects  and  seals  and  increases  the  operation  life  that  are 
vital  for  commercial  applications  [2,5].  To  speed  up  the  technology 
development,  there  have  been  a  growing  number  of  theory  and 
modeling  activities  in  this  subject  that  provide  in-depth  insights  to 
the  research  community. 

As  summarized  in  a  few  recent  papers  [6-12],  numerical  mod¬ 
els  of  varying  degrees  of  sophistication  have  been  developed. 
Broadly  speaking,  these  numerical  models  may  be  classified  as 
micro-modeling  [13-15],  macro-modeling  [8,12,16,17],  and  multi¬ 
scale  modeling  [18].  Micro-modeling  focuses  on  the  properties  of 
new  material  systems  and  membrane-electrode  assembly  (MEA) 
and  is  particularly  helpful  in  revealing  the  connection  between 
the  microstructures  and  the  properties  of  the  electrode  materi¬ 
als,  but  is  usually  incapable  of  providing  detailed  information  such 
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as  the  distributions  of  current  and  gas  species  in  a  cell  scale  and 
the  overall  power  densities  of  SOFC  stacks  in  practical  applica¬ 
tions.  Macro-modeling  aims  to  predict  the  SOFC  performance  at 
the  cell  and  stack  level  and  is  very  useful  for  improving  engineer¬ 
ing  designs,  but  is  often  limited  by  the  experimentally  determined 
model  parameters  and  does  not  pursue  the  prediction  of  the 
fuel  cell  performances  with  different  electrode  microstructures. 
Multi-scale  modeling  in  principle  combines  the  advantages  of  both 
micro-modeling  and  macro-modeling,  but  is  still  in  its  infancy  and 
integrated  three-dimensional  (3D)  models  are  yet  to  appear.  In  fact, 
most  macro-models  are  not  truly  3D  even  if  they  appear  to  be. 
For  example,  the  3D  modeling  at  the  cell  and  stack  scale  reported 
in  Refs.  [16,17]  assumed  implicitly  straight  line  paths  for  the  cur¬ 
rent  conduction  and  gas  transport  in  MEA  consisting  of  anode, 
electrolyte  and  cathode  layers  and  should  be  regarded  as  quasi- 
3D  models.  Moreover,  these  models  as  well  as  most  other  models 
ignore  issues  in  the  stack  assembly  such  as  the  effects  of  intercon¬ 
nect  ribs  on  the  performance  of  planar  SOFCs  [14,19-21  ].  Inclusion 
of  the  rib  effects  in  cell-  and  stack-level  models  is  critically  impor¬ 
tant  for  the  model  predicting  power  as  both  the  experiment  and 
theory  have  shown  that  the  performance  degradation  from  ideal 
SOFC  single  cell  to  stack  cell  is  very  significant  [19-23]. 

In  the  present  work,  multi-physics  3D  numerical  models  cou¬ 
pling  the  mass  transports  in  the  porous  electrodes,  the  electric 
current  conduction  in  the  MEA,  and  the  electrochemical  reactions 
at  the  gas-electrode-electrolyte  three-phase  boundary  (TPB)  are 
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described.  The  current  conduction  in  MEA  is  described  by  the 
partial  differential  equations  for  the  electronic  conduction  in  the 
electrodes  and  ionic  conduction  in  the  electrolyte.  The  dusty  gas 
model  is  used  to  describe  the  mass  transport  process  in  the  porous 
electrodes.  The  electrochemical  reaction  at  the  TPB  is  treated  by 
the  Butler-Volmer  equation  with  the  experimentally  determined 
parameters.  The  models  use  realistic  geometries  of  SOFC  cells  in 
practical  applications  and  consider  the  co-,  counter-  and  cross-flow 
designs  in  a  stack  assembly  as  well  as  take  the  effects  of  interconnect 
ribs  on  the  gas  transport  and  current  conduction  into  consideration. 
An  equivalent  numerical  method  that  is  mathematically  rigorous  is 
proposed  to  reduce  the  mesh  node  numbers  and  enhance  the  effi¬ 
ciency  of  numerical  simulation.  As  examples  of  applications,  the 
models  are  used  to  simulate  the  stack  cell  performance  and  to  opti¬ 
mize  the  rib  design  for  the  co-flow,  counter-flow  and  cross-flow 
SOFC  stacks. 


2.  Methods 


2.1.  Physical  model 

A  planar  SOFC  stack  is  composed  of  repeating  single  fuel  cells 
connected  in  series.  In  addition  to  the  core  part  of  a  MEA  with 
porous  anode,  dense  electrolyte  and  porous  cathode  layers,  another 
important  component  of  a  stack  cell  is  the  interconnect  with  par¬ 
allel  channels  dug  in  both  sides  to  distribute  the  gas  flow  across 
the  cell.  Examples  of  a  repeating  single  cell  with  typical  channel 
designs  are  shown  in  Fig.  1.  Close-up  views  of  the  cell  structures 
are  illustrated  in  Fig.  2.  The  typical  material,  geometric  and  oper¬ 
ational  parameters  are  shown  in  Table  1 .  As  shown  in  Table  1 ,  the 
fuel  cell  model  is  assumed  to  be  isothermal  and  all  physical  prop¬ 
erties  are  determined  at  the  temperature  of  700  °C.  Admittedly 
the  material  properties  such  as  the  electrolyte  conductivity  and 
the  activation  polarization  are  strongly  influenced  by  the  tempera¬ 
ture  field  in  a  working  SOFC  stack  [16,17].  However,  we  are  leaving 
the  model  with  coupled  thermal  conduction  equations  to  a  future 
development.  Therefore,  the  modeling  results  here  may  be  viewed 
as  the  baseline  scenario  when  the  temperature  distributions  are 
approximately  uniform.  Nevertheless,  the  following  mathematical 
descriptions  are  formulated  to  be  also  applicable  to  non-isothermal 
cells  if  coupled  with  the  thermal  conduction  equations.  Moreover, 
the  modeling  results  for  the  interconnect  rib  designs  are  expected 
to  be  valid  also  for  non-isothermal  SOFC  stacks  as  the  dominating 
physics  is  the  competition  between  the  current  collection  by  the 
interconnect  ribs  and  the  gas  transport  in  electrodes  that  are  not 
very  sensitive  to  the  temperature  field. 


Fig.  2.  Schematics  of  planar  SOFC  stack  cell  models:  (a)  co-flow;  (b)  counter-flow; 
(c)  cross-flow. 


Table  1 

Geometric,  material,  and  basic  model  parameters  for  a  planar  SOFC  cell  using  hydro¬ 
gen  fuel  and  air  oxidant. 

Parameters  Value 


Anode  thickness  (p,m) 

Cathode  thickness  (p,m) 

Electrolyte  thickness  (|jim) 

Active  cell  size  (Zx  x  ly,  cm  x  cm) 

Anode  porosity 
Anode  tortuosity 

Anode  mean  particle  diameter  (p,m) 

Cathode  porosity 
Cathode  tortuosity 

Cathode  mean  particle  diameter  (|Jim) 

Cell  temperature  (°C) 

Anode  conductivity  (s  m-1 ) 

Cathode  conductivity  (s  nrr1 ) 

Electrolyte  conductivity  (s  m-1 ) 

Diffusion  volume  of  H2  (m3  moF1 ) 

Diffusion  volume  of  H20  (m3  mol-1 ) 

Diffusion  volume  of  02  (m3  moF1) 

Diffusion  volume  of  N2  (m3  moF1 ) 

Molar  mass  of  H2  (kg  moF1 ) 

Molar  mass  of  H20  (kg  moF1 ) 

Molar  mass  of  02  (kg  moF1 ) 

Molar  mass  of  N2  (kg  moF1 ) 

Permeability  of  anode  (m2) 

Permeability  of  cathode  (m2) 

Viscosity  of  fuel  (Pa  s) 

Viscosity  of  air  (Pas) 

Knudsen  diffusion  coefficient  of  H2  (m2  s-1 ) 
Knudsen  diffusion  coefficient  of  H20  (m2  s_1 ) 
Knudsen  diffusion  coefficient  of  02  (m2  s_1 ) 
Knudsen  diffusion  coefficient  of  N2  (m2  s-1 ) 
Inlet  fuel/air  pressure,  Patm  (Pa) 

Hydrogen  molar  fraction  in  fuel 
Hydrogen  utilization 
Oxygen  utilization 

Rib-electrode  contact  resistance  (£2 cm2) 

Cell  output  voltage  (V) 


750 

50 

10 

10  x  10 
0.38 
3 
1 

0.3 

3 
1 

700 

3.356  xlO4  exp(  1392/7) 

1.223  xlO4  exp(-600 /7) 

3.34  xlO4  exp(- 10,300/7) 

7.07  x  10“6 

12.7  x  10“6 

16.6  x  10-6 

17.9  x  10“6 

2  x  10“3 

18  x 10-3 

32  x  10-3 

28  x  10“3 

7.93  x  10“16 

3.06  x  10-16 

2.8  x  10“5 

4  x  10-5 
4.37  x  10-4 
1.46  x  10“4 
7.64  x  10“5 
8.17  x  10“5 
1.013  x  105 
0.97 

70% 

20% 

0.05 

0.7 


Fig.  1.  Single  cell  unit  of  a  planar  SOFC  stack:  (a)  co-  or  counter-flow;  (b)  cross-flow. 
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The  overall  cell  performance  depends  on  the  operating  cell 
potential  and  the  output  current.  The  operating  cell  potential  (Vcell) 
can  be  formally  expressed  as  [20]: 

VceH  =  E q  —  ?7ASR;a  ~  da  ~  *7act;a  —  ?7ohm;a  —  *7ohm;el  —  ^7ohm;c 


-  lactic  -rfc-  »MSR;c  (1) 

where  Eg  is  the  Nernst  potential  when  the  partial  pressure  of  H2, 
H20  and  02  are  all  at  1  atm  (Eg  =  1.01  Vat 973 K or 700 °C),  pa ct;a 
and  lactic  are  respectively  the  anode  and  cathode  activation  over¬ 
potentials,  ?7ohm;a>  ^ohm;ei and  ^0hm;c  the  ohmic  overpotentials  in  the 
anode,  electrolyte  and  cathode,  p AsR;a  and  p Asr;c  the  anode-rib  and 
cathode-rib  interface  overpotentials  due  to  the  contact  resistance 
at  the  material  boundaries,  p°  and  p°  the  anode  and  cathode  con¬ 
centration  balance  potentials.  The  concentration  balance  potentials 
are  calculated  as  [20]: 


RT . 
ZF  n 


Ph 


2, a 


Ph2o,£ 


Itpb 


(2a) 

(2b) 


where  R  is  the  universal  gas  constant,  T  the  temperature  at  Kelvin, 
F  the  Faraday  constant,  pH2,a,  PH2o,a  and  Po2,cItpb  are  the  partial 
pressure  of  H2,  the  partial  pressure  of  H20  at  the  anode  TPB  and  the 
partial  pressure  of  02  at  the  cathode  TPB,  respectively.  The  details  of 
computing  the  other  overpotential  terms  in  Eq.  (1 )  will  be  described 
later. 

The  average  output  current  density  of  the  SOFC  cell  may  be  cal¬ 
culated  by 


I  r1*  rly 

J  lxX'y  Jx=0  iy=JZ 


dxdy 


(3) 


where  lx  (ly)  is  the  cell  length  along  the  x  (y)  direction,  jz  is  the  z 
component  of  the  current  density  flux  vector  (see  Fig.  2  for  the  axes 
definition). 


(T1-75/(p(p1//3  +  i/V3)  ))[(1/Mj)  +  (1/M,)]1/2)  the  binary  diffusion 
coefficient,  £,  r,  v,  and  Mz  the  porosity,  tortuosity,  diffusion  vol¬ 
ume  and  molar  mass  of  species  i,  respectively  [28,29].  The  required 
parameters  are  shown  in  Table  1. 

Summing  both  sides  of  eq.5a  and  eq.5b  and  using  xi  +  x2  =  1 
and  Vxi  =  -Vx2,  one  gets 


N2  Ni  Vp  /  kp  f  X\  x2  \\ 

D?r~D?rRT  v 

Substituting  Eq.  (6)  into  Eq.  (5)  yields 


Nl= _ _ VCl 

D$  +  x,D$+x2DfK 

_ D1KD2K  Xi  Vp  _  CTfcVp 

Dff  +  XiD$  +  x2Dff  RT  /i 

DeffDeff 

= _ 12  1 K _ VCl 

Deff  +  XlDeff+X2Dfff 

(  DU<D2K  k\„ 

Cl  (  RTctot  (Dg  +  X!  Dff  +  x2Df<)  +  H  J  p 


(6) 


(7) 


_  ^diffusion  _|_  jyconvective 

Therefore,  the  effective  diffusion  coefficients  of  the  species  are 


D  i  = 


DeffDeff 

U\2U\K 


^12  +  Xl  ^2I<  +  X2^1/C 

DeffDeff 

U\2U2I< 


°2  bg  +  x,D$  +  x2D?K 
The  effective  molar  flow  velocity  is  given  by 

i1 


u= ~ 


DeffDeff 

U\KU2K 


RTctot  (Dff  +  x,D^+x2DfJ) 


Vp 


(8a) 

(8b) 

(9) 


2.2.  Gas  transport  in  porous  electrode 

2.2.1.  Governing  equations 

The  mass  transport  processes  in  porous  electrodes  are  governed 
by  the  mass  diffusion  and  convective  equations.  For  species  i,  the 
transport  equation  can  be  expressed  by  [24] 

V  •  Nt  =  V  •  (— DjVCj  +  qu)  =  Rt  (4) 

where  D,,  c*,  and  Rt  are  respectively  the  diffusion  coefficient,  molar 
concentration  and  reaction  rate  of  species  i,  u  the  velocity  of  the 
mixture  flow  and  Nt  =  -Dz  vc*  +  qu  the  total  molar  flux  of  species  i. 
Di  and  u  may  be  obtained  from  the  dusty  gas  model  as  described 
below. 

2.2.2.  Dusty  gas  model 

The  original  dusty  gas  model  in  binary  gas  can  be  expressed 
[25-27] 

—iff  (p7*.  +  *1  vp + x, VP5^)  ™ 

|| +  -  ~if  (pV'2  +*lVp+*2Vp^f)  (5b> 

where  xz  is  the  molar  fraction  (xt  =  q/J2jcj^  ^  Permit¬ 
tivity,  /x  the  viscosity,  p  the  total  gas  pressure,  D^(=  sDu<l r) 
the  effective  Knudsen  diffusion  coefficients  and  D?.ff(=  sDyl r) 
the  effective  binary  diffusion  coefficients,  Dy(=  3.198  x  10-8  x 


where  the  -(/<//x)vp  term  corresponds  to  the  result  by  the  usual 
Darcy’s  law  [7,20]. 

2.3.  Electrical  conduction 

The  electronic  charge  transport  in  the  electrodes  and  the  ionic 
charge  transfer  in  the  electrolyte  are  governed  respectively  by 

-V.(ae  We)  =  0  (10a) 

— V  •  (cTj  Wj)  =  0  (10b) 

where  cre  (cr[)  is  the  electronic  (ionic)  conductivity  of  the  electrode 
(electrolyte),  Ve  (VI)  the  electric  potential  in  the  electrode  (elec¬ 
trolyte).  — cre  We(— CTj  WJ  is  the  flux  vector  of  the  electronic  (ionic) 
current  density.  The  electronic  and  ionic  potential  differences  along 
the  electrical  current  flux  paths  yield  the  ohmic  overpotentials, 
*7ohm;a»  ^?ohm;c  and  p0hm;el- 

The  electric  potential  loss  inside  the  interconnect  plate  is 
assumed  to  be  negligible  due  to  the  high  conductivity  of  the  metallic 
material.  The  local  current  densities  cross  the  interconnect/anode 
C/i^a)  and  the  cathode/interconnect  (jc-+i)  interfaces  are  deter¬ 
mined  by  the  associated  electric  potential  changes,  or  the  interface 
overpotentials: 

.  _  Ve,I/a  ~  Ve,a/I  _  PASR;a  ri1ax 

Jl^a  acd  A9R  (lla) 

xvcmcontact  xvoKContact 

,  _  Ve,c/\  ~  ^e,I/c  _  PASR;c 

Jc^\  asr  asr  (lib) 

Xvoivcontact  x\.DKContact 
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•  Experimental  data1 
■  Fitted  curve 
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0 


X. 


•i 


1 


Current  Density  (AcnT ) 

Fig.  3.  Comparison  of  the  experimental  and  theoretical  I-V  curves. 


where  Ve,i/a  and  Ve  a/I  are  respectively  the  interconnect  and  anode 
electric  potentials  at  the  anode-interconnect  boundary,  Ve,c/i  and 
VeI/c  the  cathode  and  interconnect  electric  potentials  at  the 
cathode-interconnect  boundary,  ASRcontact  the  area  specific  contact 
resistance. 


2.5.  Numerical  method 

Based  on  the  observation  that  the  coupling  between  the  gas 
channel  flows  and  the  above  described  physical  processes  is 
weak  [17],  only  the  anode,  electrolyte  and  cathode  domains  are 
included  in  the  physical  model  that  is  assumed  to  be  isother¬ 
mal  at  the  present.  The  gas  channel  flows  and  the  interconnect 
ribs  can  be  replaced  with  proper  boundary  conditions  for  gas 
transport  and  electric  conduction  [20].  For  simplicity,  the  chan¬ 
nels  (ribs)  on  the  anode-side  and  on  the  cathode-side  are  of 
the  same  width.  As  the  TPBs  are  effective  only  at  the  region 
very  close  to  the  electrode/electrolyte  interfaces  [1],  the  TPBs 
are  assumed  to  be  located  only  at  the  electrode/electrolyte 
interfaces  and  there  are  no  electrochemical  reactions  inside  the 
electrodes. 

Symmetries  are  used  to  reduce  the  model  size  and  enhance 
numerical  efficiency  whenever  possible.  For  co-  and  counter-flow 
designs,  a  multiple-channel  cell  of  the  size  of  lx  x  ly  x  lz  (x:  the  fuel 
channel  flow  direction)  may  be  represented  by  a  single  channel  por¬ 
tion  of  the  size  of  lx  x  dpitch  x  lz,  where  dpitch  is  the  pitch  width  and 
dpitch  =  ^channel  +  ^rib*  Here  dchannel  and  drib  are  one  half  of  the  chan¬ 
nel  width  and  one  half  of  the  interconnect  rib  width,  respectively. 
The  reduction  is  significant  as  ly  in  our  model  is  10  cm,  while  dpitch 
is  only  a  few  millimeters. 


2.4.  Electrochemistry  reaction 

The  current  densities  generated  by  the  electrochemical  reac¬ 
tions  at  the  anode  and  cathode  TPBs  are  described  by  the  empirical 
Butler-Volmer  equation  as  [14,27,30,31  ] 


2.5.1.  Boundary  conditions  (BCs) 

The  boundary  settings  for  the  mass  transport  equations  are 
shown  in  Table  2.  The  molar  concentrations  at  the  channel/anode 
or  channel/cathode  interface,  c°,  c°,  c°,  and  c°,  are  related  to  the 
molar  fractions  by  the  ideal  gas  equation  of  state.  The  total  gas 
pressure  at  the  electrode/channel  interface  is  set  at  1  atm.  As  the 


Jo,  a 


ja  —  < 


JO,  a 


(  Ph2,tpb 
(  Ph2,tpb 

V  PH2,f 


0.11 

0.11 


forpH2o,TPB  <  14,  000  Pa 
forpH2o,TPB  >  14,  000  Pa 


(12a) 


(12b) 


where  j0,a  and  j0iC  are  respectively  the  anode  and  cathode  exchange 
current  densities  and  may  be  determined  by  fitting  the  experimen¬ 
tal  I-Vc  urves  of  button  cells.  Fitting  the  experiment  of  Ref.  [  1  ]  for  the 
temperature  of  700 °C  gives  j0a  =  4280  Am-2  and  j0c  =  1070 Am-2 
(a  hydrogen  molar  fraction  of  0.9  in  fuel  was  used  in  the  fitting  as 
deduced  from  the  experimental  open  circuit  voltage).  As  shown  in 
Fig.  3,  the  theoretical  I-V  curve  agrees  with  the  experiment  very 
well  when  the  output  current  density  is  less  than  3  Acm-2  and  is 
sufficient  for  practical  applications.  Notice  that  the  exchange  cur¬ 
rent  densities  may  be  temperature  dependent  [32],  however,  the 
temperature  dependence  is  not  a  concern  here  as  we  are  dealing 
with  isothermal  models. 

The  activation  polarizations  in  Eq.  (12),  pact;a  and  pact;c,  are 
related  to  the  electric  and  balance  potentials  by 

Pact;a  =  ^e,a/el  “  ^i,el/a  “  Pa  (13a) 

Pact;c  =  ^i,el/c  —  ^e,c/el  —  Pc  (13b) 

where  Ve,a/ei  and  Viel/a  are  respectively  the  anode  and  elec¬ 
trolyte  electric  potentials  at  the  anode-electrolyte  boundary,  Vx  el/c 
and  \4,c/ei  the  electrolyte  and  cathode  electric  potentials  at  the 
cathode-electrolyte  boundary. 


overall  results  of  the  gas  transports  in  porous  electrodes  are  not 
very  sensitive  to  the  details  of  the  channel  flows  [17,19,20],  a  linearly 
distributed  hydrogen/oxygen  concentration  along  the  flow  is  spec¬ 
ified  on  the  channel/electrode  boundary  based  on  the  fuel/oxygen 
utilization.  For  co-flow,  both  c°  and  c°  decrease  with  x.  For  counter¬ 
flow,  c°  decreases  with  x  while  c°  increases  with  x. 

There  is  no  apparent  symmetry  in  a  cross-flow  design  that  can 
be  used  to  simplify  the  model  and  different  fuel  channels  corre¬ 
spond  to  different  c°.  Multi-physics  modeling  of  the  whole  cell 
can  be  memory  and  CPU  demanding.  Fortunately,  excessive  oxi¬ 
dant  is  used  and  the  oxygen  utilization  is  usually  low  in  practical 
applications  [17,19],  i.e.,  c°  does  not  vary  widely  over  the  whole 
cell.  The  overall  cell  performance  may  be  approximately  repre¬ 
sented  by  a  single  fuel  channel  pitch  with  the  average  c°  for  the 
boundary  condition  of  the  air  channels  (an  oxygen  molar  fraction 
of  0.19  is  used  in  this  work  and  corresponds  to  the  oxygen  uti¬ 
lization  of  20%).  Note  that  all  air  channels  and  cathode  ribs  are 
presented  in  this  simplified  model,  but  their  length  is  dpitCh  instead 

Of  ly. 

The  contact  resistance  is  set  on  the  interface  between  intercon¬ 
nect  ribs  and  the  electrodes.  The  area  specific  contact  resistances  at 
the  anode/interconnect  and  the  cathode/interconnect  boundaries 
are  assumed  to  be  uniform.  The  boundary  settings  for  the  elec¬ 
tronic  and  ionic  charge  transfer  equations  are  shown  in  Table  3. 
Due  to  the  model  difference,  the  distribution  of  electrical  bound¬ 
aries  of  the  cross-flow  design  is  also  different  from  that  of  the  co- 
or  counter-flow  design. 
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Table  2 

Boundary  settings  for  mass  transports  in  electrodes. 


Boundary 

BC  type 

Anode/channel  interface 

Anode  TPB 

All  others 

Insulation/symmetry 

(H2)  molar  concentration 

(H20)  molar  concentration 

(H2)  Inward  molar  flux 

(H20)  Inward  molar  flux 

BC 

c° 

c° 

-Ja/2/F 

Ja/2/F 

Boundary 

Cathode/channel  interface 

Cathode  TPB 

All  others 

BC  type 

(02)  molar  concentration 

(N2)  molar  concentration 

(02)  Inward  molar  flux 

(N2)  Inward  molar  flux 

Insulation/symmetry 

BC 

c° 

c° 

-Jc/4  IF 

0 

Table  3 

Boundary  settings  for  the  electronic  and  ionic  charge  transfer  equations. 


Boundary 

Rib/cathode  interface 

Rib/anode  interface 

Cathode  TPB 

Anode  TPB 

All  others 

BC  type 

Reference  potential 

Reference  potential 

Inward  current  flow 

Inward  current  flow 

Electric  insulation 

Electronic 

BC 

0  0 

uj 

1 

0 

jc 

-Ja 

Ionic 

BC 

~jc 

ja 

2.5.2.  Geometric  model  transformation 

An  anode-supported  SOFC  often  consists  of  widely  different 
geometric  dimensions.  As  indicated  in  Table  1,  the  typical  elec¬ 
trolyte  thickness  is  in  the  order  of  10  fxm,  while  the  cell  size  is 
about  10  cm  x  10  cm.  To  be  numerically  stable,  the  large  difference 
in  model  dimensions  may  require  a  mesh  with  a  great  number  of 
degrees  of  freedom.  A  reduced  dimensional  difference  is  desirable 
for  improving  the  numerical  efficiency.  Based  on  analysis  of  the 
ionic  conducting  equation,  here  we  propose  a  method  that  enlarges 
the  electrolyte  thickness  by  n  folds  with  the  corresponding  change 
in  the  electrolyte  conductivity  to  moderate  the  thin  film  effect  in 
the  meshing  step  and  decreases  the  number  of  degrees  of  freedom. 

As  described  above,  the  electrolyte  only  involves  the  partial  dif¬ 
ferential  equation  (PDE)  for  the  ionic  current  conduction  and  the 
associated  boundary  conditions.  If  the  boundary  conditions  remain 
the  same  and  the  ionic  current  is  equivalent  when  the  thickness 
of  the  electrolyte  layer  is  scaled,  the  same  physics  is  ensured.  This 
provides  the  basis  for  the  mathematical  transformation.  When  the 
thickness  of  the  electrolyte  layer  thickness  is  enlarged  by  n  folds 
while  the  x-y  dimension  remains  unchanged,  x  ->  x!  =  x,  y  ->  y'  = 
y,  z  ->  z'  =  nz,  the  same  current  flux  before  and  after  the  inhomo¬ 
geneous  scaling  requires: 

j'z=h  (i4) 

Based  on  }x  =  -ax{dV/dx),  Jy  =  -ay{dV/dy),  Jz  =  - crz(dV/dz ) 
and  V(x,y,  z)  =  V\xr ,yf ,  z'\  it  is  easy  to  know  that  the  ionic  con¬ 
ductivity  of  the  electrolyte  should  transform  with  the  geometric 
scaling  as 


In  our  3D  modeling,  the  thickness  of  electrolyte  is  scaled  from 
10  [xm  to  100  fxm,  so  the  effective  conductivity  of  the  electrolyte  is 
changed  from  isotropic  (crel,  crel,  crel)  to  anisotropic  (crel/10,  crel/10, 
lOffel)- 

The  above  mathematical  analysis  is  in  principle  general  and 
applicable  to  other  material  components  of  the  SOFC  model  and 
other  geometric  scaling  and  can  efficiently  decrease  the  mesh 
node  numbers  and  enhance  the  efficiency  of  numerical  simulation. 
Notice,  however,  the  physics  remains  correct  only  if  all  the  mathe¬ 
matical  equations  and  the  boundary  conditions  associated  with  the 
scaling  transformation  are  considered  properly. 

2.5.3.  Meshes  and  solutions 

The  finite  element  commercial  software  COMSOL 
MULTIPFISICS®  Version  3.4  [24]  is  used  in  the  present  study 
to  solve  the  PDEs  with  the  appropriate  boundary  settings.  The 


COMSOL  stationary  nonlinear  solver  uses  an  affine  invariant  form 
of  the  damped  Newton  method  [24]  to  solve  the  discretized  PDEs 
with  a  relative  convergence  tolerance  of  1  x  10-6.  Tetrahedral 
meshes  were  used  in  the  3D  models.  For  the  co-flow  and  counter¬ 
flow  models  and  with  a  pitch  width  of  2  mm  (consists  of  a  channel 
width  of  1.2  mm  and  rib  width  of  0.8  mm),  the  number  of  mesh 
elements  is  20698  and  the  number  of  degrees  of  freedom  is  96940. 
The  model  geometry  for  the  cross-flow  design  is  quite  complicated. 
For  a  pitch  width  of  2  mm,  for  example,  there  are  50  pairs  of  air 
channels  and  ribs  in  the  model  with  a  10  cm  fuel  channel.  To  reduce 
the  computational  cost,  the  10  cm  fuel  channel  model  is  subdivided 
into  five  sub-models  each  with  a  fuel  channel  length  of  2  cm.  The 
fuel  concentration  boundary  conditions  for  the  sub-models  are 
determined  based  on  the  assumed  linear  variation  description. 
The  results  for  the  five  sub-models  are  then  combined  to  produce 
the  results  for  the  whole  model.  For  a  pitch  width  of  2  mm,  the 
number  of  mesh  elements  is  44383  and  the  number  of  degrees  of 
freedom  is  209794  for  each  cross-flow  sub-model. 

3.  Results  and  discussion 

As  examples  of  applications,  the  above  described  models  are 
used  to  simulate  the  stack  cell  performance  and  optimize  the  rib 
design,  providing  the  3D  distributions  of  physical  quantities  that 
are  unavailable  with  2D  models  and  testing  the  validity  of  the  rib 
design  optimizations  by  2D  models. 

3.1.  Distributions  of  physical  quantities 

The  distributions  of  physical  quantities  such  as  the  hydrogen 
and  oxygen  molar  fractions,  current  density  and  electrical  poten¬ 
tial  in  a  stack  cell  are  important  information  for  the  evaluation  of 
a  stack  design.  The  following  results  are  obtained  with  the  basic 
model  parameters  listed  in  Table  1  and  the  stack  cell  assembly 
of  a  pitch  width  of  2  mm,  with  1.2  mm  for  the  channel  width  and 
0.8  mm  for  the  rib  width.  Due  to  the  dimensional  difference  in  the 
cell  geometry,  most  plots  use  unequal  axis  scales. 

3.1.1.  Co-  and  counter-flow  models 

Both  the  co-  and  counter-flow  stack  cells  produce  the  same 
total  current  corresponding  to  an  average  current  density  of 
0.387  A  cm-2  for  the  MEA.  It  is  interesting  to  note  that  the  aver¬ 
age  current  density  obtained  by  the  present  3D  model  is  identical 
to  that  obtained  by  a  previous  2D  model  with  the  same  model 
parameters  [20],  supporting  the  usefulness  of  a  2D  modeling  in  this 
regard.  The  3D  distributions  of  physical  quantities  for  the  co-flow 
design  are  shown  in  Fig.  4.  The  distributions  of  physical  quantities 
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Fig.  4.  The  distributions  of  physical  quantities  for  the  co-flow  cell:  (a)  hydrogen  and  oxygen  molar  fractions;  (b)  anode  TPB  current  density  by  electrochemical  reaction;  (c)  x 
component  of  the  current  density  flux;  (d)  y  component  of  the  current  density  flux;  (e)  z  component  of  the  current  density  flux;  (f)  the  electric  potential. 


for  the  counter-flow  design  are  characteristically  the  same  and  are 
not  shown. 

The  distributions  of  the  hydrogen  molar  fraction  on  the  anode 
and  oxygen  molar  fraction  on  the  cathode  are  shown  in  Fig.  4a.  As 
can  be  seen  from  Fig.  4a,  the  hydrogen  molar  fraction  at  the  anode 
TPB  is  similar  for  both  the  area  under  the  fuel  channel  and  the  area 
under  the  interconnect  rib.  The  oxygen  molar  fraction  at  the  cath¬ 
ode  TPB  is,  however,  strongly  affected  by  the  rib  presence.  Except  for 
the  region  very  close  to  the  rib-channel  boundary,  there  is  virtually 
no  oxygen  in  the  area  under  the  rib.  The  H2  and  02  distributions  for  a 
2D  (y-z)  cross-section  of  a  given  x  in  the  3D  model  and  the  underly¬ 
ing  mechanism  are  basically  the  same  as  the  previous  2D  modeling 
results  [20].  The  3D  results,  however,  have  the  definite  advantage 
of  showing  clearly  the  variations  of  H2  and  02  distributions  along 
the  flow  direction. 

Fig.  4b  shows  the  current  density  distribution  at  the  anode  TPB 
generated  by  the  electrochemical  reactions  that  is  affected  by  both 
the  H2  and  02  distributions.  Along  the  y-direction,  the  TPB  current 
density  is  dominated  by  the  02  distribution  and  there  is  basically 


no  current  generation  for  the  area  under  the  rib.  For  the  area  under 
the  channel,  the  TPB  current  density  profile  is  mainly  determined 
by  the  H2  distribution  along  the  x-direction  as  the  variation  of  the 
02  distribution  is  relatively  small. 

Fig.  4c-f  shows  the  distributions  of  current  density  flux  and  elec¬ 
trical  potential  on  the  electrodes.  The  current  density  distributions 
show  very  interesting  features.  Although  only  the  z  component  of 
the  current  density  flux  at  the  electrode-rib  interface  is  responsible 
for  producing  the  useful  output  current,  the  x  and  y  components  of 
the  current  density  flux  are  of  comparable  or  even  larger  magni¬ 
tude.  The  y  component  of  the  current  density  flux  is  necessary  for 
collecting  by  the  rib  the  current  generated  under  the  channel  area. 
The  magnitude  of  the  y  component  of  the  current  density  flux  can 
be  very  large  due  to  that  the  electrode  layers  are  thin  and  the  cross- 
section  area  for  the  y-component  flux  is  small.  The  x  component  of 
the  current  density  flux,  which  is  ignored  in  a  2D  model,  is  opposite 
to  the  fuel  flow  direction  in  the  anode  side,  but  is  along  the  fuel  flow 
direction  in  the  cathode  side.  The  magnitude  of  the  x  component  of 
the  current  density  flux  can  also  be  large  and  is  larger  than  2  A  cm-2 
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Table  4 

Average  output  current  densities  (j)  for  the  five  segments  of  the  cross-flow  models.  x°  is  the  H2  molar  fraction  at  the  channel-anode  interface  boundary. 

0.97-0.834  0.834-0.698  0.698-0.562  0.562-0.426  0.426-0.29 

]( A  cm-2)  0.482  0.421  0.376  0.337  0.298 


anode 


Fig.  5.  The  distributions  of  physical  quantities  for  the  cross-flow  cell:  (a)  hydrogen  and  oxygen  molar  fractions  at  the  TPBs;  (b)  anode  TPB  current  density  by  electrochemical 
reaction:  (c)  x  component  of  the  current  density  flux;  (d)  y  component  of  the  current  density  flux;  (e)  z  component  of  the  current  density  flux;  (f)  the  electric  potential. 
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in  the  middle  of  anode.  The  profile  of  the  current  density  flux  can  be 
explained  by  the  distribution  of  the  electric  potential  on  the  elec¬ 
trodes  shown  in  Fig.  4f.  For  example,  the  difference  of  the  electric 
potential  between  the  fuel  inlet  and  outlet  is  about  -0.01  V  and 
the  electric  conductivity  of  anode  is  about  1.4  x  105  S  m-1  at  700  °C, 
consequently  the  average  anode  current  density  along  the  fuel  flow 
direction  is  about Jx  =  A V/R  «  -0.01  V/(0.1  m/1.4  x  105  Snrr1)  = 
-1.4  A  cm-2,  agreeing  with  the  data  shown  in  Fig.  4c.  In  fact,  the 
electric  potential  distribution  is  also  affected  by  the  current  den¬ 
sity  distribution,  with  the  z  component  of  the  current  density  flux 
playing  the  dominant  role.  The  anode  electrical  potential  at  the 
anode-rib  interface  is  determined  by  the  ohmic  polarization  due 
to  the  contact  resistance  (Eq.  (11a)).  As  the  z  component  of  the  cur¬ 
rent  density  at  the  anode-rib  interface  is  higher  for  the  fuel  inlet 
than  for  the  fuel  outlet  due  to  the  hydrogen  distribution  (Fig.  4a 
and  e),  the  electric  potential  at  the  fuel  inlet  is  lower  than  that  at 
the  fuel  outlet,  causing  the  above  described  current  flow  along  the 
x  direction.  In  this  sense,  the  x  component  of  the  current  flux  is 
determined  by  the  z  component  of  the  current  flux. 

3.1.2.  Cross-flow  model 

As  described  in  Section  2.5.3,  the  whole  cross-flow  model  is  sim¬ 
ulated  with  five  sub-models  each  with  a  fuel  channel  length  of  2  cm 
and  an  average  02  molar  fraction  of  0.19  for  the  air  channel.  The 
pitch  width  is  also  2  mm,  1.2  mm  for  the  channel  width  and  0.8  mm 
for  the  rib  width.  The  obtained  average  current  densities  for  the 
five  sub-models  are  shown  in  Table  4.  The  overall  average  current 
density  for  the  cross-flow  model  is  0.383  A  cm-2,  about  1%  smaller 
than  the  result  for  the  co-  or  counter-flow  model. 

The  distributions  of  representative  physical  quantities  for  the 
cross-flow  design  are  shown  in  Fig.  5.  The  2D  distributions  of  the 
hydrogen  molar  fraction  on  the  anode  TPB  and  oxygen  molar  frac¬ 
tion  on  the  cathode  TPB  are  shown  in  Fig.  5a.  The  distribution  of 
oxygen  molar  fraction  at  the  cathode  TPB  for  the  cross-flow  design 
is  characteristically  the  same  as  that  for  the  co-flow  design,  with 
abundant  oxygen  for  area  under  the  air  channel  and  almost  no  oxy¬ 
gen  for  area  under  the  cathode  rib.  Similarly,  the  distribution  of  the 
anode  TPB  current  density  (Fig.  5b)  is  mainly  determined  by  the  dis¬ 
tribution  of  oxygen  molar  fraction  at  the  cathode  TPB  due  to  that 
the  low  ionic  conductivity  of  the  electrolyte  requires  the  oxygen 
ions  to  move  along  almost  the  shortest  possible  paths  in  the  elec¬ 
trolyte  in  order  to  minimize  the  ohmic  polarization.  This  is  the  main 
reason  for  the  very  similar  overall  current  outputs  for  the  co-  and 
cross-flow  designs  even  though  the  distributions  of  other  physical 
quantities  such  as  the  hydrogen  molar  fractions  may  appear  to  be 
very  different  for  the  co-flow  and  cross-flow  designs. 

As  shown  in  Fig.  5a,  the  hydrogen  molar  fraction  distribution 
on  the  anode  TPB  for  the  cross-flow  cell  exhibits  interesting  fea¬ 
tures  and  is  characteristically  different  from  that  for  the  co-flow 
cell.  The  main  tendency  of  the  H2  distribution  on  the  anode  TPB 
is  to  decrease  along  the  flow  direction,  but  the  profile  shows  local 
maxima  and  is  completely  different  from  the  linear  pattern  set  on 
the  channel/anode  interface.  This  is  in  fact  due  to  that  the  amount 
of  hydrogen  received  at  the  TPB  region  has  good  correspondence 
with  the  H2  distribution  in  the  fuel  channel  (similar  to  that  shown 
in  Fig.  4a),  but  the  FI2  consumption  by  the  small  current  gener¬ 
ation  at  the  anode  TPB  regions  shadowed  by  the  cathode  ribs  is 
small,  leaving  the  unconsumed  hydrogen  to  form  local  maxima  in 
the  corresponding  areas. 

Fig.  5c-f  shows  the  3D  distributions  of  the  x,  y,  and  z  compo¬ 
nents  of  the  current  density  flux  and  the  electrical  potential  on  the 
electrodes.  Like  the  H2  distribution  on  the  anode  TPB,  the  distri¬ 
butions  of  the  current  density  flux  and  the  electrical  potential  for 
the  cross-flow  may  appear  to  be  quite  different  from  that  for  the 
co-flow.  However,  the  distributions  of  the  current  flux  and  elec¬ 
trical  potential  for  the  cross-flow  may  also  be  easily  understood 


Fig.  6.  Results  of  the  rib  width  optimizations  by  3D  co-flow  models:  (a)  comparison 
of  the  optimized  rib  widths  with  the  3D  and  2D  models;  (b)  the  output  current 
densities  for  the  optimal  rib  width  and  for  drj b  =  dPitCh/2. 

like  the  H2  distribution  for  the  cross-flow.  The  mechanism  behind 
the  distributions  of  the  current  flux  and  electrical  potential  for  the 
cross-flow  cell  is  in  principle  the  same  as  that  described  above  for 
the  co-flow  cell,  namely  the  H2/02  and  TPB  current  density  distri¬ 
butions,  the  current  collections  by  the  ribs  and  the  balance  of  the 
electrical  potential,  and  is  not  repeated  here. 

3.2.  Rib  width  optimizations  in  3D  models 

A  series  of  calculations  based  on  2D  models  with  the  same  basic 
model  parameters  listed  in  Table  1  have  been  performed  for  find¬ 
ing  the  optimal  rib  widths  for  given  pitch  sizes.  The  optimal  rib 
widths  were  found  by  the  2D  models  to  be  dependent  on  the  pitch 
widths  linearly  (but  not  proportionally),  drib  =  A  +  B  x  dpitch,  and 
the  parameters  in  the  linear  relationship,  A  and  B ,  depended  prac¬ 
tically  only  on  the  contact  resistance  at  the  electrode-rib  interface 
[20].  Here  we  perform  similar  calculations  with  the  3D  models  to 
test  the  validity  of  the  2D  modeling  results. 

3.2.1.  Co-  and  counter-flow  models 

The  results  of  the  rib  width  optimizations  for  the  co-  and 
counter-flow  designs  are  practically  the  same  and  only  the  results 
for  the  co-flow  are  discussed  here.  The  optimization  results  with 
the  3D  models  for  a  few  representative  contact  resistances  are  com¬ 
pared  with  the  2D  results  and  shown  in  Fig.  6a.  As  can  be  seen,  the 
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Fig.  7.  Results  of  the  rib  width  optimizations  by  3D  cross-flow  models:  (a)  compar¬ 
ison  of  the  optimized  rib  widths  with  the  3D  and  2D  models;  (b)  the  output  current 
densities  for  the  optimal  rib  width  and  for  dr\ b  =  dPitch/2. 


results  by  the  2D  and  3D  models  are  quite  similar.  This  is  in  fact  not 
surprising  as  the  2D  model  can  be  viewed  as  a  cross  section  in  the 
3D  model  and  the  variation  of  the  H2  concentration  in  the  fuel  chan¬ 
nel  in  the  3D  model  may  be  well  approximated  with  the  average  H2 
concentration  used  in  the  2D  model  due  to  the  weak  dependence 
of  the  optimal  rib  width  on  the  H2  concentration  [20].  The  nearly 
identical  optimization  results  by  the  2D  and  3D  models  also  pro¬ 
duce  nearly  identical  optimal  output  current  densities.  To  illustrate 
the  benefit  of  the  rib  optimization,  the  optimal  current  densities 
are  compared  with  that  of  a  naive  design  of  drib  =  dpitch/2  in  Fig.  6b. 
When  the  contact  resistance  is  small  and  dpitCh  >  2  mm,  the  optimal 
output  current  densities  are  10-20%  larger  than  the  naive  design. 
When  the  contact  resistance  is  relatively  large,  the  benefit  of  the 
rib  optimization  over  the  naive  design  is  not  obvious  as  the  naive 
design  is  close  to  be  optimal. 

3.2.2.  Cross-flow  model 

The  geometry  and  the  distributions  of  physical  quantities  in 
a  cross-flow  model  are  very  different  from  that  in  a  co-flow  or 
counter-flow  model,  or  that  appeared  in  a  2D  model.  It  may  be 
rather  questionable  intuitively  for  a  2D  model  to  be  representative 
of  a  3D  cross-flow  model.  As  shown  in  Fig.  7a,  however,  the  results 
of  the  rib  width  optimizations  by  the  3D  cross-flow  model  and  the 
simplified  2D  model  are  very  close  to  each  other.  Similarly,  the  out¬ 
put  current  densities  for  the  cross-flow  cells  shown  in  Fig.  7b  are 
also  close  to  that  for  the  co-flow  cells  shown  in  Fig.  6b,  though  the 


Fig.  8.  Comparison  of  the  theoretical  and  experimental  I-V  results.  The  difference 
between  the  dashed  line  and  the  dash-dotted  line  shows  a  possible  performance 
improvement  of  the  stack  cell  if  a  better  rib-channel  structure  is  used. 


former  may  be  1  -2%  smaller  than  the  latter.  The  ability  of  a  sim¬ 
plified  2D  model  to  represent  a  3D  cross-flow  model  is  in  fact  not 
surprising  due  to  the  dominance  of  the  effects  of  the  cathode  ribs 
on  the  TPB  current  generation  and  the  inconsequential  role  of  the 
anode  ribs  on  the  hydrogen  transport  and  the  TPB  current  genera¬ 
tion,  as  described  above  in  Section  3.1.2.  The  results  are  gratifying 
since  the  optimization  modeling  is  important  for  improving  the 
engineering  design  and  the  2D  modeling  is  much  more  economic 
than  the  3D  modeling. 

3.2.3.  Comparison  with  the  experimental  data 

An  example  has  been  given  in  Fig.  3  to  show  that  the  theoret¬ 
ical  model  may  reproduce  the  experimental  I-V  data  by  adjusting 
the  anode  and  cathode  exchange  current  densities  that  are  depen¬ 
dent  on  a  variety  of  microstructural  parameters  [15].  The  theoretical 
model  is  further  used  here  to  explain  the  different  I-V  results  of 
a  single  cell  and  an  identical  cell  in  a  three-cell  stack  observed 
experimentally  [22].  As  shown  in  Ref.  [14],  the  current  collec¬ 
tor  of  platinum  mesh  used  in  the  single  cell  testing  [22]  may  be 
described  by  model  parameters  of  dpitCh  =  0-3  mm,  drib  =  0.05  mm 
and  ASRcontact  =  0.008  C2  cm2.  Fig.  8  shows  the  comparison  of  the 
experimental  data  of  Ref.  [22]  for  the  temperature  of  700  °C  and  the 
theoretical  results  withj0a  =400  Am-2  and  j0c  =  100  Am-2  (based 
on  the  experimental  open  circuit  voltage  of  1.2  V,  a  hydrogen  molar 
fraction  of  0.996  in  the  input  fuel  was  used  in  the  fitting,  correlating 
well  with  the  pure  hydrogen  fuel  used  in  the  experiment).  As  shown 
in  Fig.  8,  the  agreement  between  the  theoretical  and  experimental 
I-V  curves  is  satisfactory. 

The  experimental  three-cell  stack  consists  of  five  rib-channel 
pitches  in  a  square  cell  of  5  cm  x  5  cm  [22],  corresponding  to 
dpitch  =  5  mm  and  drib  =  0.7  mm  in  the  theoretical  model.  The  input 
fuel  flow  rate  may  be  deduced  from  the  fuel  utilization  for  a  given 
output  current  density  stated  in  the  experiment  [22].  Based  on 
these  data,  the  theoretical  performance  of  a  single  cell  in  the  three¬ 
cell  stack  may  be  obtained  and  shown  together  with  experimental 
results  in  Fig.  8.  Notice  that  the  experimental  current  densities  in 
Ref.  [22]  were  scaled  by  a  factor  of  18.49/25  here  as  the  fuel  cell  area 
used  in  the  experimental  presentation  was  18.49  cm2  instead  of  the 
actual  25  cm2  used  in  this  work.  Clearly,  the  experimental  observa¬ 
tion  of  the  decreased  performance  of  a  stack  cell  in  comparison 
with  the  corresponding  single  cell  is  quantitatively  explained  by 
the  theory,  demonstrating  the  validity  and  the  predictive  power  of 
the  theoretical  model. 
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Notice  that  the  pitch  width  of  5  mm  used  in  the  experiment 
is  unnecessarily  large  and  detrimental  to  the  stack  cell  perfor¬ 
mance  [19].  If  a  pitch  width  of  2  mm  and  a  rib  width  of  0.42  mm 
were  used  in  the  stack  construction,  the  stack  cell  performance 
would  be  improved  notably,  as  shown  in  Fig.  8.  For  example,  the 
output  current  density  of  such  a  stack  cell  for  the  operating  volt¬ 
age  of  0.7  V  is  expected  to  be  0.37  A  cm-2,  an  increase  of  about 
30%  over  the  result  of  0.29  A  cm-2  for  the  experimental  stack 
cell. 

4.  Summary 

We  have  presented  three-dimensional  multi-physics  numerical 
models  for  planar  SOFC  cells  with  the  co-,  counter-  and  cross-flow 
stack  designs.  The  models  are  capable  of  handling  the  electronic 
conduction  in  the  electrodes,  ionic  conduction  in  the  electrolyte, 
mass  transport  in  the  porous  electrodes  and  electrochemical  reac¬ 
tions  on  the  three  phase  boundaries  as  well  as  the  critical  role 
of  interconnect  ribs  on  the  cell  performance.  A  geometric  scal¬ 
ing  algorithm  that  is  in  principle  general  is  proposed  to  improve 
the  meshing  efficiency  of  the  3D  model.  Numerical  examples  for 
different  flow  designs  are  presented  and  the  3D  distributions  of 
physical  quantities  are  displayed  and  analyzed  with  interesting  fea¬ 
tures  that  are  unavailable  to  the  2D  models.  The  models  are  also 
used  to  find  the  optimal  rib  widths  for  given  channel-rib  pitch 
widths.  The  results  by  the  3D  models  show  that  the  stack  cell 
performances  and  the  optimal  rib  widths  are  all  very  similar  for 
co-,  counter-,  and  cross-flow  designs  and  can  be  well  represented 
by  the  corresponding  2D  models,  simplifying  the  future  modeling 
work  as  far  as  the  cell  current  output  and  the  rib  optimization  are 
concerned. 
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